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Abstract 

The numerical properties of algorithms for finding the intersection of 
sets depend to some extent on the regularity of the sets, but even more 
importantly on the regularity of the intersection. The alternating pro- 
jection algorithm of von Neumann has been shown to converge locally 
at a linear rate dependent on the regularity modulus of the intersection. 
In many applications, however, the sets in question come from inexact 
measurements that are matched to idealized models. It is unlikely that 
any such problems in applications will enjoy metrically regular intersec- 
tion, let alone set intersection. We explore a regularization strategy that 
generates an intersection with the desired regularity properties. The reg- 
ularization, however, can lead to a significant increase in computational 
complexity. In a further refinement, we investigate and prove linear con- 
vergence of an approximate alternating projection algorithm. The anal- 
ysis provides a regularization strategy that fits naturally with many ill- 
posed inverse problems, and a mathematically sound stopping criterion 
for extrapolated, approximate algorithms. The theory is demonstrated 
on the phase retrieval problem with experimental data. The conventional 
early termination applied in practice to unregularized, consistent prob- 
lems in diffraction imaging can be justified fully in the framework of this 
analysis providing, for the first time, proof of convergence of alternating 
approximate projections for finite dimensional, consistent phase retrieval 
problems. 

1 Introduction 

The role of local regularity for nonconvex minimization problems or nonmono- 
tone variational inequalities is well-established. In broad terms, a generalized 
equation is said to be "regular" (or "metrically regular" ) if the distance from a 
proposed solution to an exact solution can be bounded by a constant multiple 
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of the model error of the proposed solution. A particular focus has been the 
proximal point algorithm and alternating projections |Tir5 HTTl[TS] ■ 

It is often the case, however, that the problems in question are ill-posed] 
in other words, there is no constant of proportionality between the model error 
and the distance of an approximate solution to the true solution. For some algo- 
rithms such an ill-posedness would not prevent the iterates from converging to 
a best approximate solution, but numerical performance will suffer. An example 
of such behavior can be observed with the classical alternating projection algo- 
rithm of von Neumann [22] applied to a general feasibility problem: that is, the 
problem of finding the intersection of sets. Ill-posedness for feasibility problems 
can be characterized by problem inconsistency, that is, the nonexistence of an 
intersection of the sets in question. More generally, the feasibility problem will 
be ill-posed if the intersection vanishes under arbitrarily small perturbations of 
the sets. 

For the applications we have in mind, at least one of the sets in question 
comes from a finite precision measurement or calculation. It is quite reasonable 
to expect an inconsistency between the idealized model and the measured data, 
which can be represented as a perturbation of the idealized data set. When only 
two convex sets are involved, alternating projections can be shown to converge 
to nearest points [B] Theorem 4] , however the rate of convergence will in general 
be arbitrarily slow. For other algorithms ill-posedness leads to instability in the 
sense that the iterates do not converge to a fixed point. The Douglas Rachford 
algorithm, for example, applied to inconsistent feasibility problems has no fixed 

points [1H2Q3]. 

Insofar as ill-posed problems can be regularized, the theory cited above can 
be applied to numerical methods for the regularized problems. Our focus here 
is on a particular regularization for ill-posed feasibility problems and efficient 
approximate projection algorithms. The problem of nonconvex best approxima- 
tion was considered in 03E] where the focus was on instability of the Douglas 
Rachford algorithm resulting from problem inconsistency. A relaxation of this 
algorithm was proposed that has fixed points for inconsistent problems and has 
been successful in practice [15]. As is often the case for relaxed projection al- 
gorithms, there is no systematic rule for choosing the relaxation parameter. It 
was shown in |14) that the size of the relaxation parameter at the solution is 
related to the optimal gap distance between the sets. This observation suggests 
a different approach to algorithmic design that is based on regularization of the 
underlying problem rather than stabilization of the algorithm as was the focus 
in Q3]. 

We further develop this viewpoint here, where we study local regularization 
of the underlying problem while retaining the character of the original problem. 
In particular, we expand one of the sets in order to create an intersection with 
all the desired regularity properties described in [IT]. The strategy is a local 
regularization in the sense that indicator functions are still used as the central 
penalty function, in contrast to |14) where the indicator function was relaxed to 
a distance function. One then can apply any number of algorithms for finding 
the intersection of regularized sets. We are particularly interested in projection 
algorithms and specifically the classical alternating projection algorithm. Wc 
show in section [3] that, for the problems of interest to us, such a regularization 
of the sets results in a significant increase in the complexity of computing the 
corresponding projections. To address computational complexity of the regu- 



2 



larized problem we consider approximate alternating projections based on the 
projection operators of the original, unregularized problem. An approximate 
algorithm is stated in sectional We prove local linear convergence of this algo- 
rithm to a solution of the regularized problem under regularity assumptions that 
are natural for regularized problems. In section[5]we apply a specific approxima- 
tion motivated in section [3] to the approximate projection algorithm and prove 
that this approximation is guaranteed to succeed under certain conditions. We 
demonstrate the effectiveness of this approximation in section[S]with an example 
from diffraction imaging with real experimental data. We do not claim that the 
approximate alternating projection algorithm is the best, or even a very good 
strategy for solving this particular problem. However to our knowledge, our 
analysis yields the first mathematically sound stopping criteria for alternating 
projections applied to the phase retrieval problem. Our goal is to demonstrate 
the theory and to motivate the adaptation of our proposed regularization and 
approximation to more sophisticated projection algorithms. 

2 Notation, Definitions and Basic Theory 

We begin with basic theory and notation. For the most part, we present only 
the results with pointers to the literature for interested readers. The setting we 
consider is finite dimensional Euclidean space E. The closed unit ball centered 
at x is denoted by B(x); when it is centered at the origin, we simply write B. 
We denote the open interval from a to 6 by (a, b) ; the closed interval is denoted 
as usual by [a, b}. 

Given a set C C E, we define the distance function and (multivalued) pro- 
jection for C by 

dc(x) = d(x,C) = inf{||z - x\\ : z G C} 
Pc{%) — argmin {\\z — x\\ : z G C}. 

If C is closed, then the projection is nonempty. Following 17, Definition 1.6] 
we define the normal cone to a closed set C G E as follows: 

Definition 2.1 (normal cone) A vector v is normal to a closed set C C E at 
x, written v G Nc(x) if there are sequences x k — > x and v — > v with 

v k G {t(x k - z) 1 1 > 0, z € P c (x k ) } for all keN. 

The vectors v k are proximal normals to C at z G Pc{x k ) and the cone of 
proximal normals at z is denoted Nq{z). 

It follows immediately from the definition that the normal cone is a closed 
multifunction: for any sequence of points x k — > x in C, any limit of a sequence 
of normals v k G Nc(x k ) must lie in Nc(x). The relation of the projection to 
the normal cone is also evident from the definition: 

zePc(x) => x-zeNc(z). (2.1) 

Notice too that N c (x) = {0} x e int C. 
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Definition 2.2 (basic set intersection qualification) A family of closed sets 
C\,C<z^... C m C E satisfies the basic set intersection qualification at a point 
x € HjCi, if the only solution to 

m 

$Jfc=°> Vi&N Ci {x) (i = 1,2,..., m) 

8=1 

is yi = /or i = l,2,...,m. PFe saj/ that the intersection is strongly regular at 
x if the basic set constraint qualification is satisfied there. 

In the case m = 2, this condition can be written 

iv Cl (s)n-iv C2 (5) = {0}. 

The two set case is is called the basic constraint qualification for sets in |17[ Def- 
inition 3.2] and has its origins in the the generalized property of nonseparabil- 
ity [T|5] which is the n-set case. It was later recovered as a dual characterization 
of what is called strong regularity of the intersection in [101 Proposition 2] . This 
property was called linear regularity in [11] . The case of two sets also yields the 
following simple quantitative characterization of strong regularity. 

Proposition 2.3 (Theorem 5.16 of [llj ) Suppose thatC\ andC2 are closed 
subsets o/E. The intersection C\ l~l Ci satisfies the basic set intersection quali- 
fication at x if and only if the constant 

c := max{(u, v) | u G N Cl (x) n B, t)6-JV Cj (l)nl}<l. (2.2) 



Definition 2.4 (angle of regular intersections) We say that the intersec- 
tion C\ n C2 is strongly regular at x with angle 9 := cos _1 (c) > where c is 
given by (I2.2p . 

In order to achieve linear rates of convergence of alternating projections to 
the intersection of sets, we require pointwise strong regularity of the intersection 
[llj . In the absence of this property the above definitions suggest a general 
regularization philosophy: promote strong regularity. This is most obviously 
achieved by augmenting at least one of the sets by some e ball: C\{e) = C\ + 
eB, for instance. Similar ideas been used extensively in the development of 
proximally smooth sets by Clarke, Stern and Wolenski 0. We pursue this idea 
in section [3] with the generalization that the ball, or "tube" around the set of 
interest is with respect to a generic distance in the image space of a continuous 
mapping, the tube having no relation to the native space in which the projectors 
onto the sets are defined. 

Somewhat stronger results are possible when the sets have additional regu- 
larity. We call a set C C E is prox-regular at a point x G C if the projection 
mapping Pq is single- valued around x [19) . Convex sets, in particular, are 
prox-regular. More generally, any set defined by C 2 equations and inequalities 
is prox-regular at any point satisfying the Mangasarian-Fromovitz constraint 
qualification, for instance. 
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Proposition 2.5 (angle of normals of prox-regular set) Suppose the set 
C C E is prox-regular at the point x G C . Then for any constant 5 > 0, any 
points y,z € C near x and any normal vector v £ Nc(y) satisfy the inequality 

(v, z-y) < 5\\v\\ ■ \\z - y\\. 

Proof. This is a special case of the same property for super regular sets ( 
Definition 4.3] and [TTJ Proposition 4.4]) since by (TTJ Proposition 4.9] prox- 
regularity implies super regularity. 

Alternatively, for prox-regular sets we can proceed directly from [191 Propo- 
sition 1.2] which shows that, for any sequences of points y k ,z k g C converging 
to x and any corresponding sequence of normal vectors v k £ Nc(y k ), there exist 
constants e, p > such that 



-v , z 



' 1 y k )<l\\z k -y k t 



for all large k. Since for any fixed 5 > we will eventually have \\z k — y k \\ < —, 
it follows that 

<«*, z k - y k ) < S\\v k \\ ■ \\z k - y k \\ 

for k large enough. □ 

The next result builds upon Proposition l2.5l and provides bounds on the angle 
between sets in the neighborhood of a point in a strongly regular intersection 
of a closed and a prox-regular set. In [TTJ Theorem 5.2] implications (|2.3p and 
(|2.4j) are used to characterize sets for which linear convergence of the alternating 
projections algorithm holds. We do not seek such generality here and are content 
with identifying classes of sets which satisfy these conditions, namely prox- 
regular sets. The proof of the following assertion can be found in the proof of 
Theorem 5.16 of [TTj . 

Proposition 2.6 Let M, C C E fee closed. Suppose that C is prox-regular at a 
point x € MflC and that M and C have strongly regular intersection at x with 
angle 8. Define c := cos(#) and fix the constant d with c < d < 1. There exists 
a constant e > such that 

x e M n (x + eM), ue-N M {x)nm\ , 
yeCn{x + eM), veN c {y)r\B j [u,v)sc, [z.6) 



and, for some constant S £ [0, —^-) 
y,z € C n (x + eB] 

v € N c (y)nn 



(v, z - y) < 6\\z - y\\. (2.4) 



In what follows, we define an approximate alternating projection algorithm 
in terms of the distance of the normal cone associated with the approximate pro- 
jection to the "true" normal cone. In order to guarantee that for our proposed 
approximation we can get arbitrarily close to the true projection, we need the 
notion of convergence of the associated normal cone mappings. Let S : E Y 
denote a set- valued mapping where Y is another Euclidean space. We define 
the domain of S to be the set of points whose image is not empty, that is 

domS':={x \ S(x) =^0}. 
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Following Definition 4.1] wc define continuous set- valued mappings rela- 
tive to some subset D as those which are both outer and inner semicontinuous 
relative to D. 

Definition 2.7 (continuity of set-valued mappings) A set-valued mapping 
S : E =4 Y is continuous at a point x G D relative to D a K if 

S(x) C 

{y | V x k 5> x, 3 A" > suc/i t/iaf /or k > K, y k -> y wzt/i y fc G S^) } 

and 

{y | 3 x fe ^> 3 y fc -> y wit/i y fc G S(x*) } C S(z) 

where ^> indicates that the sequence lies within D. We denote this as S{x) — > 
S(x) for all sequences x ^> x. 

3 The problem 

In this section we formulate our abstract problem and motivate the regulariza- 
tion and approximation strategies that we propose. Our initial, naive problem 
formulation involves finding points x G C C E, a Euclidean space, that explain 
some measurement b G Y modeled as the image of the continuous mapping 
g : E -> Y , that is 

Find x eC(l M 

for 

M := {x G E |ff(x) = 6}. 

The set C usually captures a qualitative feature of solutions, such as non- 
negativity, or a prescribed support. If 6 is a physical/empirical measurement, 
it is likely that the intersection is empty, or that the solution consists only of 
extremal points. In the case of measurements with discrepancies modeled by 
statistical noise, the noise could be Gaussian or Poisson distributed (among 
still other possibilities) . To accommodate a variety of instances we consider the 
following regularizations of the set Mo: 

M e := {x G E | d^gix), b) < e} (3.1) 

where e > and d^ is a Bregman distance defined by 

d<t>{z,y) ■= <t>{z) - <t>{y) - 4>'(y)(z - y) 

for cj> : Y — > (—oo, +oo] strictly convex and diffcrentiablc on int (dom (f>) . The 
Bregman distance with <j> := ^|| • || 2 corresponds to the Euclidean norm which 
is appropriate for Gaussian noise. If Y = R m and 

m ( tlogt -t for t > 

<p(y) = KVj) for h W : = { for t = 

i =1 \+oo fort<0 

then the Bregman distance leads to the Kullback-Leibler divergence, 

m 

d<t>{z, y) = KL(x, y) := ^ z j lo S ~ + Vj ~ z r ( 3 - 2 ) 
3=1 Vl 

The Kullback-Leibler divergence is appropriate for Poisson noise. 
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Remark 3.1 The regularization p. II) bears some resemblance to closed neigh- 
borhoods of the type X{e) :— {x \ d(x,X) < e} considered by Clarke, Stern and 
Wolenski [7] in their development of proximally smooth sets, except that the 
neighborhood around the set of interest is with respect to a generic distance in 
the image space of a continuous mapping, the neighborhood having no relation 
to the metric upon which the projectors onto the sets are defined. Still, we 
will rely on prox-regularity of the regularized set for the approximation strategy 
discussed in Section [5J □ 

Regardless of the distance, the first algorithm we consider for finding this 
intersection is the classical alternating projection algorithm. 

Algorithm 3.2 (exact alternating projections) Choose x° <E C . For k = 

1, 2, 3, . . . generate the sequence {x 2k } C C with x 2k G Pcix 2 ^ 1 ) where 
the sequence {x 2k+1 } consists of points x 2k+1 <G PM c ( x2k )- 

We show next that the projection onto the fattened set M e could be considerably 
more costly to calculate than for the unregularized set Mq. This motivates the 
approximate projection algorithm studied in section 2] 
We want to compute 

x* € Pm c (x) '■= ar g mm :reAf e \\\ x " x \\ 2 - 

Assume d<f,(g(x), b) > e, then we seek a solution on the e-sphere around b with 
respect to d^. This is an instance of a trust region problem. 

Suppose that x G P^ e (x) and that the standard constraint qualification 
holds, that is 

- Vd (.g(x,6))*77 = O, V >0 =^ 77 = 0. (3.3) 

Then 

(S-S) + Vd (g(S),&)*ry = (jj > 0) (3.4) 
d^g(x),b)-e = 0. (3.5) 

These are the standard KKT conditions (see, for example [101 Theorem 10.6]). 
Numerical methods for computing the projection Pu t (x) involve solving a pos- 
sibly large-scale nonlinear system of equations with respect to x and 77; this 
could well be as difficult to solve as the original problem. 

Example 3.3 (affine subspaces) Let E = W\ Y = R m with m < n. Take 
g to be the linear mapping A : R" —> R m and d ( j ) (x,y) — ^\\x — y\\ 2 (that 
is, (j)(x) = 7}\\x\\ 2 ). The projection can then be written as the solution to a 
quadratically constrained quadratic program: 

minimize nlla; — z\\ 2 
subject to \\\Ax - b\\ 2 < e. 

For small problem sizes this can be efficiently solved via interior point methods. 
Still, even the most efficient numerical methods cannot compare to computing 
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the projection onto the afhne space Mq := {i£E \Ax = b} which has the 
trivial closed form 

P Mo (z) = (/- A T {AA T )- 1 A)z + A T (AA T )~ 1 b. 

This suggests an alternative strategy for computing the projection onto the 
"fattened" set. 

Indeed, we can efficiently compute the projection Pm, (z) as a convex com- 
bination of the points y — Pm {%) and z 

x* = \ e z + (1 — X e )y 

where A e e [0, 1) solves ^(1 — A) 2 ||z — y\\ 2 = e. This also has a closed form: 
the quadratic formula. For general Bregman distances such shortcuts are not 
available, but this forms the basis for our approximations. □ 

Example 3.4 (boxes) Let E = Y = R". Define g : R" -> R" by 

g(x) = (\x!\ 2 , \x n \ 2 ) T 

and, again, let the distance d$ be the standard normalized squared Euclidean 
distance to some point b € K™. The projection can then be written as the 
solution to the nonconvex optimization problem 

minimize ills — x\\ 2 

subject to \ X^id^jl 2 — ^i) 2 = e - 

Notice that the corresponding set M e is not convex: the origin is projected in 
the positive and negative direction in each component. Generally, nonconvex 
problems are hard to solve. On the other hand, the projection onto the box 
with length 2b, y = (j/i, . . . , y n ) € Pm (x), is trivial and has the form 

See [5] for analysis of this projection in higher dimensional product spaces. 

For this example there is no shortcut to computing the projection Pm { for 
e > 0, but we show below that the convex combination of the projection of a; onto 
Mq and x is a effective approximation that still yields linear rates of convergence 
for the method of alternating projections for finding the intersection of M e n C. 
□ 

Remark 3.5 We note that in both of the above examples the constraint qualifi- 
cation p.3[) is no longer satisfied in the limit e = for the set M e . This obviously 
does not prevent us from calculating the projection onto the set Mq. Indeed, as 
we showed, the projection sometimes even has an explicit representation. □ 
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4 Inexact alternating projections 



There is more than one way to formulate inexact algorithms. One template 
for this is to add summable error terms to the operators involved in the exact 
algorithm. Another approach - the one we take here - is less general but has 
a more geometric appeal. More to the point, it is appropriate for our intended 
application. 

Given two iterates x 2k ~ Y £ M and x 2k £ C, a necessary condition for the 
new iterate x 2k+1 to be an exact projection on M, that is x 2k+1 £ PM(x 2k ), is 

H^fc+l _ x 2 fe|| < ||x 2fe _ ^fc-lji x 2k _ x 2k+l g NM{x ^+ly 

In a modification of (TT| we assume only that we choose the odd iterates x 2k+1 
to satisfy a relaxed version of this condition, where we replace the second part 
by the assumption that the distance of the normalized direction of the current 
step to the normal cone to M at the intersection of the boundary of M with 
the line segment between x 2k+1 and x 2k is small. 

Consider the following inexact alternating projection iteration for finding 
the intersection of two sets M, C C E. 

Algorithm 4.1 (inexact alternating projections) Fix 7 > and choose 
x° £ C and x 1 £ M . For k = 1, 2, 3, . . . generate the sequence {x 2k } C C 
with x 2k £ Pc(x ~ 1 ) where the sequence {x 2k+1 } C M satisfies 

||z 2fe+1 -z 2fe || < \\x 2k - x 2 ^ 1 ]], (4.1a) 
x 2k+1 =x 2k i£xl k+1 =x 2k , (4.1b) 



and d 

for 



N M {^)(n<i ( 4 - lc ) 



X 

and 



z k 



2k+l _ p fn~ 2k \ 
— "MD{x 2k -rz k , t>0}\ x ) 




if x 2k+1 7^ x 2k 



Note that the odd iterates x 2k+1 can lie on the interior of M. This is the major 
difference between Algorithm 14.11 and the one specified in [TT] where all of the 
iterates are assumed to lie on the boundary of M. We include this feature to 
allow for extrapolated iterates in the case where M has interior. Extrapolation, 
or over relaxation, is a common technique for accelerating algorithms, though 
its basis is rather heuristic. Empirical experience reported in the literature 
shows that extrapolation can be quite effective (see [3 [21]). The algorithm 
given in Theorem 15.11 below explicitly includes extrapolation. Our numerical 
results at the end of this paper do not contradict the conventional experience 
with extrapolation. 

Lemma T4.2I and Theorem 14.41 below were sketched in [TH Theorem 6.1] for 
the variation of Algorithm 14. ll iust described. 



Lemma 4.2 Let M, C C E be closed. Suppose that C is prox-regular at a point 
x £ M n C and that M and C have strongly regular intersection at x with angle 
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9. Define c :— cos(#) and fix the constants c with c < c < 1 and 7 < v 1 — c 2 . 
T/ien i/iere is an e > smc/i £/ia£ t/ie iterates of Algorithm ^, i] satisfy 

•'' JA ' T " 2 1 => || x 2fc + 2 - x 2fc+1 || < ryUx 2 ^ 1 - z 2fe || (4.2) 



i^+i _ <| 



/or 77 = c-^/l — 7 2 + 7VI — c 2 < 1. 

Proof. Fix c' with c < c' < c < 1 and define S — ^ (77 — 77') where 77' = 
c' \/l — 7 2 + 7\/l — c' 2 . (5 > since, as is easily verified, cyl — 7 2 + 7\/l — c 2 
increases monotonically with respect to c on [0, 1].) Since C is prox-regular at 
2; and the intersection is strongly regular, by Proposition 12.61 for this S there is 
an e > such that implications (|2.3I) and (|2.4I) hold. We apply this result here. 
The assumptions and the triangle inequality yield 

\\x 2k -x\\ < ||a; 2fe -.T 2fc + 1 || + ||x- a; 2fc+1 || <e. (4.3) 

By the definition of x 2k+1 we have x 2k+1 = (1— A)o; 2fc +A:r 2fc+1 for some A G [0, 1] 
so that 

\\4 k+1 ~x\\ = \\\{x 2k+1 -x) + (l-\)(x 2k -x)\\ 

< A||:r 2fc+1 -z|| + (l-A)||a: 2fe -x|| 

< A^ + (1-A)e<e (AG [0,1]) (4.4) 



where the last inequality combines the left hand side of (|4.2[) and (|4.3[) . Next, 
by the triangle inequality and the definition of the projection 

\\x 2k+2 -x\\ < \\x 2k + 2 -x 2k+1 \\ + \\x-x 2k+1 \\ 

< ||x 2fc -:c 2fe + 1 || + ||5:- a ; 2fc+1 || < e. (4.5) 

If x 2k+1 = x 2k then this is a fixed point of the algorithm and the result is trivial. 
Similarly, if x 2k+1 — x 2k+2 : then x 2k+1 G Cfl M and by the first condition in 
(|4.1j) this is a fixed point of the algorithm. So we assume that x 2k+1 7^ x 2k and 
define w G N M {x 2k+1 ) with \\w\\ = 1 and u := ^ll+lzlTk+i\\ ■ Now applying 
Proposition [2~B1 to x 2k+1 satisfying (|4.4|l with -w G -N M {x 2k+1 ) n B and to 
a; 2fc+2 sa ti s fyi n g (|4.5[) with — u G A^c (x 2fc+2 ) n B we have that for e small enough 

(w, u) = (—w, —u}< c . (4-6) 

In other words the angular separation between the unit vectors w and u is 
bounded below by arccosc'. 
On the other hand, define 

x 2k _ x 2k+l 



Our goal is to obtain a lower bound the angle between z and u. If it were the 
case that x 2k+1 G PAf(x 2k ) then z = w and c' would already be our bound. 
But since x 2k+1 only approximates the projection, we must work a little harder. 
Since the iterates satisfy (|4.1[) . for some w G NM(x 2k+1 ) we have \\w — ?|| < 7. 
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There are two cases to consider. If z = 0, then we are done. Otherwise z 
has length one, and 

||Hl 2 + l-7 2 
2|H 

Maximizing the left hand side as a function of ||iu|| € [1 — 7, 1 + 7] yields the 
largest possible angular separation from z, that is 

(w, z) > Vl-7 2 (4-7) 

where w = 

Note that 7 < \f\ - c 2 < Vl - c' 2 for c' < c so that c' < ^/l - 7 2 . Thus, 
combining (14.61) and (|4.7[) . we have 

(w, z) > \J\ — 7 2 > c' > (ui, u) 

arccos (i?, z) < arccos(\/l — 7 2 ) < arccosc' < arccos (u>, u) . 
It follows immediately, then, that 

< arccosc' — arccos(\/l 7 2 ) 

< arccos (w, u) — arccos (w, z) < arccos (u, z) 

which is equivalent to 

(z, u) < cos ^arccosc' — arccos(-\/l — 7 2 )^ = c'\/l — 7 2 + 7V 1 — c' 2 < 1. 

Letting 7/ — c' — 7 2 + 7V 1 — c' 2 and removing the normalization yields 
^ a .2fc_ a .afc+i j ^H-2_ s 2^i^<^|| a; 2fc_ :i .afc+i|||| a .2fc + 2_ :i .afc+i|| < (48) 

Now by our choice of e, implication (I2.4[) holds for x 2k and x 2fe+2 eCfl {x + eB} 
with — u e iV c (a; 2fe+2 ) n B, namely 

x 2fe - z 2fe + 2 ) < <5||x 2fe + 2 - z 2fe || 

which is equivalent to 

/„2fc+2 J2.k+l 2k+2 2k\ ^ rll 2/c+2 2fc||| 2fe+2 2fe+l|| 
\X — X , X — X / < 0\\X — X \\\\X — X || . 

By the triangle inequality and the dchnition of the projection 

U^fe+2 _ < u^fc+2 _ ^fc+ljl + U^fc+l _ < gj^fc+l _ x 2fe|| 

so that 

(x 2fc + 2 - z 2 ^ 1 , a; 2fe + 2 - x 2fc ) < 2<5|| 2 ; 2fc + 1 - a; 2fe || \\x 2k+2 - x 2k+1 \\. (4.9) 

Adding (gU) and gj]) yields 

H^^-x 2 ^ 1 !! 2 < {28 + r 1 l )\\x 2k+1 ~x 2k \\\\x 2k+2 -x 2k+ \ 

which by our construction of 5 yields 

|| a .afc+2_ a; 2fc + l|| < v \\ x ™+i_ x 2k^ 

as claimed. □ 
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Lemma 4.3 With the same assumptions as Lemma choose x° and x 1 so 
that 

\W -x\\ < \\x°-x\\ =/3< ^-^e (4.10) 



where e is chosen to satisfy (|4.2p . Let rj = cy 1 — 7 2 + 7VI c 2 . T/ien /or oZZ 
fc > 

||x 2fc+1 -x|| < 2^ 1 ~^ +1 <|, (4.11a) 

H^fe+i _ ^feji < ^fe < £ and (4.11b) 

|| a .2fc+2_ a .M+l|| < ^fc-K (4nc) 



// in addition M is prox-regular at x, then for all k > 

_ „fe+i 

- < 

1-7? 2 



x k+1 -x\\ < 2/3 LjllL < £ ; (4.12a) 



||x fe+1 -x fe || < /3r; fe < I and (4.12b) 
l^fe+2 _ a .fe+i|| < ^fc+i. ( 4 _ 12 c) 

Proof. The proof is by induction. For the case A: = inequality (|4.11al) holds 
trivially. Inequality (I4.11b[) follows from the triangle inequality and (|4.10p . 
Inequality (I4.11c|) then follows from (I4.11a[) . (|4.11bl) and Lemma |4~21 Since for 
the case k = inequalities (|4.11a[) - (l4.11cj) are equivalent to (|4.12aj) - (l4.12cj) this 
case is true whether M is prox-regular or not. 

To show that these relations hold for k + 1 with AI not prox-regular, note 
that by gT]) ||x 2fc+3 -2; 2 ' c + 2 || < \\x 2k+2 - x 2k+1 \\. In light of (|4TTc)) this implies 

\\x 2k+3 - x 2k+2 \\ < (3rj k+1 < ^. (4.13) 

This together with (I4.11aj) and (|4.11cj) . yields 

_-|| < - 2 ; 2fc + 2 || + || a; 2fe + 2 - a ; 2 ' £ + 1 || + ||x 2fc + 1 -x|| 

< 07i k + 1 +j3r) k+1 +2l3 1 ~ r,k+1 

1 -77 

1 — n k+1 1 — ?i fe + 2 f 

< 2/3 v k+1 +2l3—^ = 2/3—^ <i. (4.14) 

1 — 77 1 — 77 2 

Now, Lemma [O] applied to (14.13)) and (|4.14p yields 

|| x 2 fc +4 _ x 2k+3^ < _ x 2fc+2|| < ^fe+2^ (4 lg) 

As (14.13p - (|4.15p are just (|4.11a[l - (|4.11c[) with fc replaced by k+ 1, this completes 
the induction and the proof for the case where M is not prox-regular. 
If we assume, in addition, that M is prox-regular, then by (I4.12cp 

\\x k+2 - x k+1 \\ < pr, k+1 < |. (4.16) 
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This together with (|4.12aj) yields 

| + ||z fe+1 -z|| 

< (3n k+1 +[3 



\x k+2 - x\\ < \\x k+2 - x k+1 \\ + \\x k+1 -x\\ 



1 - n k+1 
l — rj 



1 — ri k+2 7) 

= t^hr^l (4 - 17) 

Now Lemma l4~2l with the rolls of C and M reversed, together with (I4.12c[) yields 

Again, since (I4.16p - (I4.18I) are just (|4.12al) - (|4.12cl) with k replaced by k + 1, this 
completes the induction and the proof. □ 

Theorem 4.4 (convergence of inexact alternating projections) Let M, C C 

E and suppose C is prox-regular at a point x G M fl C. Suppose furthermore 
that M and C have strongly regular intersection at x with angle 6. Define 
c := cos(6*) < 1 and fix the constants c € (c, 1) and 7 < %/ 1 c 2 . For x° and x 1 
close enough to x, the iterates in Algorithm ^. 1\ converge to a point in M fl C 
with R-linear rate 



\J c\Jl - 72 + j\/T^^ < 1. 



//, in addition, M is prox-regular at x, then the iterates converge with R-linear 

rate 

cy/l - 7 2 + 7\/l - c 2 < 1. 

Proof. We prove in detail the case where M is not assumed to be prox-regular. 
Choose x° and x 1 so that (|4.10p holds with e is chosen as in Lemma [4.21 Let 
n = c\/\ — 7 2 + 7VI — c 2 . To establish convergence of the sequence we check 
that the iterates form a Cauchy sequence. To see this, note that for any integer 
k = 0, 1, 2, . . . and any integer j > 2k, by (|4.11b[) and (14. 1 lc|) of Lemma |4~3"1 we 
have 

\w-x 2k \\ < ^rii^+wii 

i=2k 

< (3(r/ k + 2n k+1 + 2n k+2 + . . . ) 



1-T) 



Similarly, it can be shown that 



_ x 2fc+l|| < p f]_ 



k+1 



1-r, 



So the sequence is a Cauchy sequence and converges to some ieE. The fixed 
point of the sequence must belong to M n C and satisfies 



0„ / a 1 + r l 



\\x-X U \\ <P 
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Moreover, for all j = 0, 1, 2, 



\\x-x j \\ < /V /2 ^r|- 

We conclude that convergence is R-linear with rate ^/fj as claimed. 

The proof for the case where M is also prox-regular at x proceeds analogously 
using inequalities (|4.12al) - (|4.12c[) of Lemma 14.21 instead. □ 

Note that the worse the approximation to the projection, the slower the 
convergence. As we showed in the previous section, the projection onto the 
unfattened set can be easier (sometimes much easier) to compute than the pro- 
jection onto the fattened set, so although the rate of convergence suffers from 
taking only an approximate projection, we gain in the per-iteration complexity 
of calculating the projections. 



5 Approximate alternating projections onto fat- 
tened sets 

Theorem 5.1 Let E and Y be Euclidean spaces, and (f> : Y — > (— oo,+oo] be 
Isc, strictly convex and differentiable on int (dom <p) . Let C C E be closed and 
M t (1C^0 for all e > where M e is defined by 

M e :={xeE\f :=d4g{x),b)<e} 

for dcf, a Bregman distance to b G Y and g : E-)Y continuous with range(g) C 
dom <fi and 

l imi nf^#^>0. 

z|->oo \x\ 

Suppose that there is a 9q > and, for all e > small enough, a point xq G 
MoflC with nearby points x e at which the intersection M e (lC is strongly regular 
with angle 9 e > 8q > 0. Suppose further that C and M e are prox-regular on a 
neighborhood ofxo and that M e has nonzero proximal normals at all boundary 
points within this neighborhood. Define cq := cos(#o) < 1 and fix the constants 
c £ (co, 1) and 7 < y/l — c 2 . Compute the sequence {x k } via A Igorithm \4-l\ with 
the odd iterates generated by 

x 2k+1 = (1 - X k )x 2k + \ k xf +1 (5.1) 

for £Q fc+1 € PM (x 2k ) and X k > 0. For x° and x 1 close enough to Xq, there 
exist {Xk} > and e > such that for all k G N the iterates satisfy (|4.1aj) . and 
(I4.1cl) with {x 2k+1 } G M t , and the sequence of points converges to a point in 
M e n C with at least R-linear rate 

Cy/l - 7 2 + 7V 1 - C 2 < 1. 

The odd iterates of the proposed algorithm do not necessarily lie on the surface 
of the regularized set M e , but could be on the interior of this set. Were we 
computing true projections, all the odd iterates would lie on the boundary 
of M e - instead we take larger steps than the projections would indicate. In 
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this sense, the algorithm defined in Theorem 15.11 is a regularized approximate 
alternating projection with extrapolation. The theorem does not tell us what 
such extrapolation buys us, but at least it says that we will not do any worse 
than without it. 

We begin next developing the groundwork for the proof of Theorem 15.11 

Lemma 5.2 (level-boundedness) Let E and Y be Euclidean spaces, and <p : 
Y — > (-co, +oo] be lsc, strictly convex and differ entiable on int (dom <p). Define 
the function f := d,p(g(-),b) where d<p(y,b) is the Bregman distance of y to the 
point b £ dom <f> and the function g : E — > Y is continuous with range(g) C 
dom <fi and satisfies 

liminf ^ >0. (5.2) 

\x\— >oc \x\ 

Then the lower level sets o//, {i£E f( x ) < a } for fixed ael, are compact. 
In particular, the set argmin / is nonempty and compact and inf / = min/ > 0. 

Proof. For easy reference we recall the definition of the Bregman distance: 

f(x) := d*(jg(x), b) = 4>(g{x)) - 4>{b) - (</>'(&), g(x) - b) . 

Since range(g) C dom <f> and b G dom 4> there is an i e E at which f(x) < oo. 
Moreover, since <j) is convex, the Bregman distance is bounded below by 0, hence 
inf / > and / is proper (that is, not everywhere equal to infinity, and does 
not take the value -co on E). Also / is lsc as the composition of the sum of 
a lsc function <f> and a linear function ((f)' (b), ■) with a continuous function g. 
The lower level sets of / are therefore closed (see for instance [UJ Theorem 
1.6]). The coercivity condition (15. 2[) then implies that the lower level-sets are 
bounded [3U1 Corollary 3.27], thus the lower level sets are compact and argmin / 
is nonempty and compact. □ 



Theorem 5.3 (continuity of the level set mapping) Let f := d^,(g(-),b) 
with <f>, g, b and d^ as in Lemma \5.S\ The corresponding level-set mapping 

M(a) := {x e E \f(x) < a} (5.3) 

is continuous on [e, oo) where e := min/. 

Proof. By Lemma T5.2I M (•) is compact and dom M(-) — [e, oo) C [0, oo). Con- 
sequently the graph of M(-) is closed (in fact, closed-valued) in E x K and 
satisfies 

{y | 3 a k -> a, 3 y k -> y with y k e M{a k ) } C M(a) for all 5eR. (5.4) 

On the other hand, the inverse of the level-set mapping (the epigraphical profile 
mapping) 

M~ l {x) := {a e K \a > f{x)} 

maps open sets to open sets relative to [e, oo), that is M~ l (0) is open relative 
to [e, oo) for every open set O C E. Thus by [501 Theorem 5.7] the level set 
mapping satisfies 

M(a) c (5.5) 
(y V a k 57, 3 K > such that for k> K, y k -> y with y k e M (a k ) | 
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for all a > e. Since the right hand side of (|5.5I) is a subset of the left hand side 
of (|5.4p we have equality of these limiting procedures, and thus continuity of 
M(-) on [e, oo) according to Definition 12.71 □ 

Proposition 5.4 Let f := dj,(g(-),b) with <f>, g, b and be as in Lemma 
\5.!H\ and let M(a) be defined by f|5 . 3[) . For {a k } C [e, oo) with a k — > a where 
Z := min/. the corresponding sequence of projections onto M(a k ), PM(a k )> 
converges graphically to Pm(u), that is 



Proof. Since M(a k ) — > M(a) by Theorem 15.31 graphical convergence of the 
projection mapping follows from a minor extension of [201 Proposition 4.9] (see 
[201 Example 5.35]). □ 

In light of the discussion in section [21 our numerical strategy for approx- 
imating the projection to the regularized set M e defined by (|3.1I) will be to 
compute the intersection of the boundary of M e with line segment between the 
current iterate and the projection onto the unregularized set Mq. Specifically, 
for x $l M e we define xq = Pm (x) and calculate the point 

x e :— (l — T e )x + T e xo where r £ := min{r > | (1 — t)x + txq £ M e }. (5.6) 

The next proposition shows that this approximation can achieve any specified 
accuracy for sets with a certain regularity. This will then be used to guaran- 
tee that the approximation to the projection given by (|5.6[) satisfies (|4.1c[) on 
neighborhoods of a fixed point of Algorithm 14. II 

Proposition 5.5 (uniform normal cone approximation) Let e > and 

M e (e € [0, e]) be defined by (37Q). Let x a G M , and (x + pB) n (E \ Mi) ^ 
for p > fixed. In addition to the assumptions of Lemma \5. 6 A suppose that M e 
(e € [0, e]) is prox-regular at all points x € (xq + pB) flM e with nonzero proximal 
normals at points x € [(xq + pB) fl M e ] \ int (M e ). Then given any 7 > there 
exists an e' £ (0,e] such that for all e £ (0,e'] 



holds where zq — Pm (z), z e is given by (|5.6p and z is any point near (xq + 
pM)DM e . 

Proof. Since for all e g [0,e] the sets M e are prox-regular on (xq + pM) n 
M e , all nonzero proximal normals to M e can be realized by an r-ball on open 
neighborhoods of points on (xo + pE) n M e for r small enough [TH Theorem 
1.3.f]. There is thus a ball with radius r E > on which the nonzero proximal 
normals to M e can be realized uniformly on (xo + pM) n M e . Also by assumption, 
the proximal normal cones to all points on the boundary of (xq + pM) n M e are 
nonzero. Thus, by Definition 12.11 the normal cone to M e at all points on the 
boundary of (xq + pB) n M e can be identified with the projection of points z in 
a r e -neighborhood of this boundary. The result then follows from Proposition 
15.41 identifying the level set mapping M(e) with the parameterized set M e . □ 



gphP M(Qfc) -> gphP M (a). 




(5.7) 
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Remark 5.6 We conjecture that the assumption of prox-regularity and non- 
triviality of the proximal normal can be relaxed. The assumptions of Lemma 
15.21 are used to guarantee graphical convergence of the projection mappings; 
the issue here is that the points on the boundary of the M e generated by (|5.6p 
do not have to correspond to projections. Prox-regularity, and more restric- 
tive still, the nontriviality of the proximal normals to M e on the boundary is 
used, in essence, locally to guarantee the reverse implication of (|2.ip . Definition 
12.11 only relies on the existence of sequences of proximal normals whose limits 
constitute the normal cone. Our approximation scheme (|5.6[) . in contrast, gen- 
erates a specific sequence of points, which could conceivably correspond only 
to zero proximal normals without further assumptions on the regularity of M e , 
though we are unaware of a counterexample. That prox-regularity alone is not 
enough to assure that the proximal normal cone is nonzero is nicely illustrated 
by the set M = |a; € R 2 X2 > x^ 5 j which is prox- regular at the origin, but 

has only a zero proximal normal cone there (see [201 Fig- 6-12.]). Obviously, 
such regularity will depend on the distance d$ and the mapping g used in the 
construction of M f . □ 



Proof of Theorem \5.1\ We show first that there are > such that for any 
e > the iterates x 2k+1 lie in M c and satisfy (I4.1a[) . Consider Afc = 1 for all k. 
Then x 2k+1 = x 2 a k+1 6 M e for all k and all e > and by the definition of the 
projection 

|| x 2fc+l_ < ^ x 2k_ x 2k-l^ 

which suffices to prove the claim. 

For existence of e > such that (|4.1cp is satisfied, note that for all e suffi- 
ciently small c e := cos(6> c ) < cos($o) -—cq<1 so that the choice of c £ (cq, 1) 
satisfies c € (c £ , 1) and consequently a fixed 7 < \fl — c 2 suffices for all e suf- 
ficiently small. The result then follows immediately from Proposition 15.51 The 
assumptions of Theorem 14.41 then apply to guarantee linear convergence, which 
completes the proof. □ 



Remark 5.7 The theorem above guarantees convergence of Algorithm @TT] with 
approximation strategy given by (|5.ip for instances where the intersection of the 
unregularized problem need not be strongly regular. When the unregularized 
problem is inconsistent the strategy may fail. In particular, suppose that Mq n 
C = 0. Then for some e the intersection M e n C = for all e < e. If 7 is such 
that (|4. lc|) is only satisfied for e < e, then the proposed approximation will fail. 

To the degree that the coupling between the regularization parameter e and 
7 is weak, we can still obtain positive results. One instance where the coupling 
is very weak is if the fattened set has interior and x is some point in this 
interior. In this case c = in (|2.2I) . 7 can be arbitrarily close to 1 and the 
condition (|4.1c[) is almost trivial to satisfy. This is indeed the case for our 
intended application. Of course, the closer 7 is to 1, that is, the worse our 
approximation of the true projection, the slower the convergence; so the trade 
off between efficient computations and rates of convergence must be balanced. 
The addition of extrapolation to the approximate algorithm is meant to mitigate 
any adverse effects of the approximation. The effectiveness of extrapolation is 
illustrated in the following section. □ 
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Figure 1: Diffraction image of real object. 

6 An example from diffraction imaging 

We present an application of the theory developed here to image reconstruction 
from laser diffraction experiments produced at the Institute for X-Ray Physics 
at the University of Gottingen. Shown in Figure [T] is the observed diffraction 
image produced by an object resembling a coffee cup that has been placed in 
the path of a helium-neon laser. The imaging model is 

\Fx\ 2 = b (6.1) 

where b £ M™ is the observed image intensity, F is a discrete Fourier transform, 
| • | 2 is the componentwise (pixelwise) modulus-squared, and x £ C n is the object 
to be found. The image is corrupted by noise modeled by a Poisson distribution. 
In the context of (|3.ip the solution we seek lies in the fattened set 

M £ := {x \KL{\Fx\ 2 ,b) < e} (6.2) 

for KL(x,y) the Kullback-Leibler divergence given by (|3.2[) . This set can be 
shown to be prox-regular everywhere with nonzero proximal normals at all 
points on the boundary. To this, we add the qualitative constraint that the 
object is nonnegative (that is, real) and lies within a specified support: for a 
given index set J C {1,2, ... ,n} 

C :={x£ K™ \ Xj =0 for j £ J}. 

This set is not only prox-regular, but in fact convex. 

Despite the good features of these sets, the problem is inconsistent /ill-posed. 
The set C is a set of real vectors, but the observation b is corrupted by noise. 
If b is not symmetric, as happens to be the case here, then the image cannot 
come from a real- valued object. Sometimes practitioners will "preprocess" the 
data by symmetrizing the raw data. If this is done, then the corresponding 
feasibility problem is provably consistent, and the results of Theorem 15.11 can 
be applied. In the numerical examples below, however, we choose to keep closer 
to the true nature of the experiment and demonstrate the success of Algorithm 
14.11 as prescribed by Theorem 14.41 despite the absence of guarantees that the 
condition (|4.1c[) is satisfied. 

The state of the art for iterative methods for solving this problem can be 
found in [15] . The main problem for these algorithms is the absence of a stop- 
ping criterion. Often what is done in practice is one algorithm (often the Dou- 
glas Rachford algorithm or variants [2"ll3"ll 13]) is used to get close to a solution, 
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iteration 



Figure 2: Reconstruction and behavior of odd and even iterates of unregularized 
(data set M given by (|6.2[1 ) exact alternating projections applied to the diffrac- 
tion imaging problem. Only 500 iterations are shown. The algorithm appeared 
to find a best approximation pair after about 24,000 iterations. (x 2k —> x 2k+2 
but \\x 2k — £ 2fe+1 || is bounded above zero.) 



and then alternating projections is used to refine the image according to the 
"eye-ball" norm. In the application literature alternating projections is often 
known as the "Error Reduction" algorithm. Different communities have differ- 
ent opinions as to what constitutes a stopping criteria, but in our reading of the 
application literature, seldom do the proposed criteria involve iterates approach- 
ing a numerical fixed point. Typical behavior of alternating projections onto 
the unregularized problem, together with the corresponding reconstruction are 
shown in Figure [21 The true object was a coffee cup, which can be seen, upside 
down, in the lower right hand corner of the reconstruction in Figure [2l with 
the handle on the left hand side. The reconstruction of the true object is only 
unique up to rotations, shifts and reflections. This is why the reconstruction is 
upside down relative to the true object. 

Next we apply Algorithm 14.11 with the approximate projection computed as 
in Theorem l5.1l for different regularization parameters e and different step-length 
strategies. Figure [3] shows the reconstruction and behavior of iterates for e = 1.9 
and Afe chosen so that the iterates remain on the surface of the M e set. Figure 
2] shows the reconstruction and behavior of iterates for e — 1.2 and A& = 1 
for all k. Figure [5] shows the apparent convergence rates for different values 
of the relaxation parameter e in (|6.2j) and different settings for the step-length 
parameters A&. The black solid line shows again the change between the even 
iterates of the unregularized, exact alternating projection algorithm. The blue 
and green dashed lines show the apparent rate of convergence of the regularized 
problems without extrapolation, that is, A& is computed so that the iterates lie 
on the surface of the set M e (to numerical precision). As expected, the lower 
the value of e, the poorer the (asymptotic) rate of convergence since the sets 
are closer to ill-posedness for smaller regularization values. The red dashed- 
dotted line shows what can be gained by extrapolation. Here the step-length 
parameter A^ = 1 for all k and the algorithm proceeds with a convergence rate 
indicated by Theorem 14.41 but then terminates finitely as it finds a point on 
intersection interior to the regularized set M e . Note that the only difference be- 
tween this implementation and the unregularized exact alternating projections 
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Figure 3: Reconstruction and behavior of odd and even iterates of regularized 
(data set M e given by (|6.2p with e = 1.9) inexact alternating projections with 
Afe chosen so that the iterates lie on the surface of the M e set. 




Figure 4: Reconstruction and behavior of odd and even iterates of regularized 
(data set M e given by f|6.2[) with e = 1.2) inexact extrapolated alternating 
projections with = 1 for all k. The algorithm terminates at the 5th iterate 
which achieves condition (|4.1bl) to numerical precision. 
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Figure 5: Comparison of implementations of Algorithm I4.ll with the approx- 
imate projection computed as in Theorem I5.ll for different parameters e and 
step- length strategies (Afc) for the fattened set M e given by (|6.2I) . The black 
line is the unregularized alternating projection algorithm with exact projections. 
The blue and green lines are the regularized approximate alternating projection 
algorithms with step lengths Afc computed so that the iterates lie on the sur- 
face of the M e set. The red line is the extrapolated approximate alternating 
projection algorithm with = 1 for all k. 

implementation (the black solid line) is early termination of the algorithm. This 
is what is usually done heuristically in practice. What this example shows is 
a mathematically sound explanation of this practice in terms of regularization, 
extrapolation and approximate alternating projections. 

For this example it is not possible to compute an a priori rate of convergence 
as specified by Theorem 14.41 since the set M e has no analytic form and we are 
unable to compute the angle of the intersection. We observe a linear convergence 
rate, at least to the limit of machine precision. To a certain extent, this is beside 
the point. The value of the theory outlined above lies not with the computation 
of rates of convergence, but rather with the provision of regularization strategies 
and corresponding stopping rules. We can, however, verify numerically whether 
a point lies in the interior of the regularized set at the point of intersection 
with the qualitative constraint set. For the extrapolated example shown in 
Figure [4] it was verified that this point lies in the interior of the set M e by 
perturbing the point slightly and verifying that it still lies in the set M e . Thus, 
even if the unregularized problem is not consistent as required by Theorem 15.11 
to guarantee that the approximate projection achieves a sufficient accuracy for 
linear convergence, since the fixed point of the algorithm is an interior point, 
the required accuracy for the approximate projection is quite easy to satisfy as 
discussed in Remark lixJl 

Finally, note that while the rate of convergence for the more regularized 
problems is better, as indicated by comparing the reconstructions in Figures [3] 
and 21 the reconstruction can be poorer since this reconstruction is apparently 
further away from the ideal solution than the less regularized reconstructions. 
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7 Conclusion 



The main achievement of this note is not our algorithm. Indeed, the regularized 
extrapolated (A& = 1 for all k) inexact projection algorithm specified in Theo- 
rem l5.1l in fact has been used successfully for decades in diffraction imaging with 
heuristic stopping criteria and early termination effectively serving as the regu- 
larization. What the analysis here provides, for the first time, is a regularization 
strategy that fits naturally with many ill-posed inverse problems, and a mathe- 
matically sound stopping criterion. The conventional early termination applied 
in practice to the unregularized problem can be justified fully in the framework 
of this regularization strategy together with approximate projections. While all 
of the regularity assumptions on the sets M c and C are satisfied for the finite 
dimensional phase problem discussed in Section [51 since the unregularized phase 
problem with noise is still inconsistent, Theorem 15.11 does not apply. If exact 
projections onto the regularized sets were computed, then Theorem 5.16 of [TT] 
would suffice to prove convergence of exact alternating projections applied to 
the regularized problem. Proof of convergence of the inexact algorithm with 
extrapolation strategy Afe = 1 for all k for the regularized, inconsistent phase 
retrieval problem (what has in fact been applied in the application literature 
for decades) hinges on verifying that condition 14. lcl of Algorithm 14. II is satisfied 
locally for all iterates. This is an open problem. 
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